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ABSTRACT 

Gas has been detected in a number of debris disk systems. This gas may have arisen from grain 
sublimation or grain photodesorption. It interacts with the surrounding dust grains through a number 
of charge and heat exchanges. Studying the chemical composition and physical state of this gas can 
therefore reveal much about the dust component in these debris disks. We have produced a new code, 
ONTARIO, to address gas emission from dusty gas-poor disks around A-F stars. This code computes 
the gas ionization and thermal balance self-consistently, with particular care taken of heating/cooling 
mechanisms. Line emission spectra are then produced for each species (up to zinc) by statistical 
equilibrium calculations of the atomic/ionic energy levels. For parameters that resemble the observed 
P Pictoris gas disk, we find that the gas is primarily heated by photoelectric emission from dust grains, 
and primarily cooled through the C II 157.7 /xm line emission. The gas can be heated to a temperature 
that is warmer than that of the dust and may in some cases reach temperature for thermal escape. 
The dominant cooling line, C II 157.7 //m, should be detectable by Herschel , while the O I 63.2 /^m 
line will be too faint. We also study the dependence of the cooling line fluxes on a variety of disk 
parameters, in light of the much improved sensitivity to thermal line emission in the mid/far infrared 
and at sub-millimeter wavelengths provided by, in particular, Herschel, SOFIA, and ALMA. These 
new instruments will yield much new information about dusty debris disks. 
Subject headings: Stars: Circumstellar Matter, Stars: Planetary Systems: Formation, Stars: Planetary 

Systems: Protoplanetary Disks, Scattering, Stars: Individual: Constellation Name: 

/3 Pictoris 



1. INTRODUCTION 

During the T Tauri phase, models have shown that gi- 
ant planet formation rely heavily on the tim escale and 
mann er in which gas disks disperse (see, e.g. Ilda fc Lml 
l2004f ). We now know that some of this gas is incorporated 
into gaseous giant planets, while some is lost to the stars, 
and that most of t he gas disappears on timescales of a 
few Myr (e.g. jZ uckcrman ct al. 1995; Hai sch et al]|2001l: 
iJavawardhana et al.. .2006; .Pascucci et al.ll2006[ ). How- 
ever, we still do not know whether most of the gas disks 
around T Tauri stars is photocvaporated by the star 
(see, e.g.JG orti & HoUcnbach 2009), blown away by wind 
(|Lovelace e t al. 20 081), or accreted onto the st ar due to 
viscous diffusion ( Lvnden-Bell fc Pringlel[l974[ ). Which 
mechanism dominates will determine how long gas re- 
sides in different parts of the disk. Observing protoplan- 
etary disks near their final stages may yield clue to this 
puzzle. 

Debris disks are circumstellar disks that show emis- 
sion from d ust. These disks are detectable for up to 
a few Gyr (IHolland et all [19981 ISoangler et all [2OOII: 
iRieke et al.l 120051 )^ Both the level of emission and the 



fraction of stars with a detectable excess decay with 
time and can be interpreted, at least for the case of 
A stars, as steady state coUisional processing of plan- 
etesimals formed during the proto-pla,netary disk phase 
pominik fc DecinI [20031: iWvatt et al.) [20071 ). There is 
some evidence that gas may persist in debris disks for an 
extended period of time, at least in some systems. For 
example, atomic gas in Keplerian rotation has been de- 
tected around the ^ 12 Myr old main-sequ ence star /3 Pic 
(jHobbs et al.lll985l: iBrandck er et ai]|2004D Gas has also 
been found around post-T Tauri stars (e.g. lRedfield et all 
2007fl. and even s tars as old as 200 Myr (e.g. crHer; 
Chen &: Jural [2003h . However, it is not clear whether 
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this gas is primordial (left-over proto-planetary material) 
or secondary, generate d by, e.g. , sublimation of plan- 
etesimals ( B_eust fc ValironI I2007D. va porization of col- 
liding planetesimals (|Czechowski fc M ann 2007), and/or 
photo-desorption of dust (jChen et al.. ,2007.) . If gas is 
indeed secondary (as dynamical ar guments suggest, see 
[Fernandez. Brandeker. fc Wull200 o) , it presents an excit- 
ing new tool to study the compositions of solid bodies in 
extra-solar systems, unaccessible otherwise. In addition, 
even a small quantity of gas may also significantly affect 
the dynamics of the dust and cause grain concentrations, 
mimicking the brigh t rings seen in young debris di sks like 
HR4796A (IKlahr fc Lin 200,4 IBeslafcWi[[2007[) . Fur- 



thermore, gas present during late stages of planet for- 
mation could damp eccentricities of the planetesimals 
enough to reduce their relative speed, thereby aiding the 
build up of rocky planets. 

Unfortunately, the information currently available on 
the manner in which gas disperses in circumstellar disks 
is scarce, owing to the great difficulty in detecting gas in 
these environments. Accretion traces only the sub-AU 
region of the disk, leaving the outer gas disk largely un- 
constrained. Molecular emission from CO in young cir- 
cumstellar disks has been detected in the ne ar-infrared 
(NIR) and in the mm (e.g. lNajita et al.ll2003( ). The NIR 
emission is limited to hot CO (~ 1 000 K), while the de- 
tection of mm emission from cold CO is ha mpered by lim- 
ited sensitivity due to beam dilution (e.g. iPascucci et al.l 
12006). Direct detection of H2 in the ultra violet (UV) 
and NIR is also made possible by cit her very hot or fluo- 
rescent gas (pumped by Lya, Herczeg et al.112006 '). The 
rotational lines H2 S(0) and S(l) emitted from colder 
gas (< 100 K) are unfortunately weak and difficult to 
observe: a Spitzer IRS survey of young stars that are ex- 
pected to have gas disks yielded only upper lim its on the 
gas disk mass (upper limits > 100 M^. lChen et al. 2006). 
An alternative is to observe stellar light scattered reso- 
nantly by the gas, but this requires a spatially resolved 
disk to be detectable. To date, only t he disk around 
/3Pic has been detected in such a way (jOlofsson et al.l 
l2001l :rBrandekcr et al. 2004). Should the disk be oriented 
edge-on, the circ umstellar gas c an also be observed in 
absorption (e.g., iRedfieldl l2007t iRoberge fc WeinbergeiH 
12001 . 

Here, we study infrared atomic cooling lines as promis- 
ing agents for characterizing gas in tenuous disks, in- 
cluding late stage protoplanetary disks and debris disks. 
Three of the most abundant species, oxygen, carbon, and 
silicon, all have ground-state fine transitions at infrared 
wavelengths that act as effective coolants. Oxygen, in 
the low-density environment of these disks, is likely to 
be in the atomic form: the molecular repositories for 
oxygen (H2O, OH, and CO; e.g.:Kamp.& Bcrtoldi 2000) 
are dissociated easily in the circumstellar environment, 
while the ionization potential for oxygen is sufficiently 
high for it to remain neutral. Moreover, neutral oxygen 
does not have strong resonant transitions in the optical 
and UV, and is thus not easily rer noved by radiation pres - 
sure, unlike some other species (iFernandez et al.l 120061 ). 
Of all thermally excited cooling lines, the O I 63 /im 
line is expected to be the brightest for a large range of 
disk m asses in late-type stars (G-K: [Gorti fc Hollenbacbl 
l2004f ). However, in optically thin disks around early- type 
stars (F or earlier), carbon is expected to be significantly 
ionized (jFernandez et al. 2006), boosting the C II 158 /im 
cooling line. Finally, if cooling by the O I and C II lines 
saturates (i.e. the line fiuxes no longer rise with tempera- 
ture due to the saturated occupation of the excited level, 
see eq. [5]), Si Il-cooling may become important. 

The O I 63 /xm and C II 158 /im lines both lie within 
the spectral range of the far-infrared to sub-millimeter 
Herschel space telescope, scheduled to start science oper- 
ations in early 2010. The telescope has a passively cooled 
3.5 m mirror and three instruments PACS, SPIRE, and 
HIFI that together cover the 55-672 /im with both imag- 
ing and spectroscopy capabilities. As such, Herschel is 
well positioned to investigate the tenuous circumstellar 



disks we are discussing here and partly motivated this 
study. 

Predicting the expected luminosity of gas and inter- 
preting potential detections require detailed modeling of 
the physical circumstances, such as the radiation field 
from the star, the spatial distribution and content of the 
gas and dust, etc. Debris disks are known to be dusty 
with little gas, meaning that the set of assumptions and 
approximations that can be used are different from other 
environments, such as the interstellar medium where 
well-tested codes exist (e.g. cloudy, Ferland 2005). De- 
bris disks are generally optically thin, which simplifies 
the radiative transfer and enables modeling of more gen- 
eral three-dimensional distributions with less computa- 
tional resources. On the other hand, a more detailed 
treatment is required for the thermal balance, as fiuores- 
cence is generally important due to the strong radiation 
field from the star in combination with low gas densities. 
This implies computing the statistical equilibrium (SE) 
for a large number of energy levels in cooling species. 

We have developed a numerical model Ontario ( "op- 
tically thin thermal and ionization equilibrium code" ) to 
investigate the importance of O I 63 /im and C II 158 /im 
fine-structure cooling lines as tracers for gas in dusty 
but optically thin disks around early-type (A-F) stars. 
This is intended to comple ment similar analysi s mad e 
for later-type stars (G-K) bv lGorti fc HoUenbachl (|2004D . 
The code takes stellar flux, gas and dust profiles as input 
parameters and performs self-consistent computation of 
ionization and thermal balance (§[2]) in the disk. Once 
elemental ionizations and gas temperature have been de- 
termined, a full SE computation is performed for a num- 
ber of atomic species to identify the dominant emission 
and absorption lines. 

To test the code, we use the well-studied /3 Pictoris 
disk as a benchmark case in §[3l and compare the numer- 
ical model with observed gas disk properties. In iJ lA.U 
we keep stellar parameters at their /3Pic values but re- 
place the gas and dust distributions with more simpli- 
fied Gaussian profiles, and take this configuration as the 
fiducial debris disk. Fiducial case temperature and ion- 
ization profiles are analyzed in S I4.1.1I and ii l4.1.3l We 
also compute integrated luminosities for the main fine- 
structure cooling lines fJi l4.1.2p and identify the dominant 
emission and absorption lines fii l4.1.4[) . In §[42] we inves- 
tigate the dependence of emitted line luminosities on the 
model parameters, in particular how they scale with disk 
gas, geometry, mass and composition, and the spectral 
type of the central star. Conclusions are drawn in §[51 

Throughout the paper, we adopt the cgs units except 
where noted. 

2. THERMAL PROCESSES MODELED IN ONTARIO 

The ONTARIO model was designed to simulate gas 
in dusty disks around early-type stars. The code self- 
consistently computes the thermal and ionization states 
of the gas in this environment. All gas is assumed to 
be present in atomic or ionized form: the hard radi- 
ation fields surrounding these early-type stars rapidly 
photo-dissoci ate molecules into individual atoms or ions 
(jKamp et al]|2003l : iJonkheid et al.ll2006[ ): thus no chem- 
istry computation is included. Ionizations for atomic el- 
ements from hydrogen to zinc are considered, up to the 
second ionization state. The model assumes both dust 



and gas to be optically thin, so no complex radiative 
transfer calculations are performed. This makes the code 
effectively zero-dimensional as each bin is treated inde- 
pendently of all others, greatly increasing computational 
performance. Details of the code, including how ioniza- 
tion and statistical equilibrium is treated, is presented in 
the Appendix. Here, we discuss the main processes that 
enter the thermal balance. 

A number of heating and cooling mechanisms that de- 
termine the gas temperature in debris disks have been 
investigated by Kamp fc van Za dclhofi (2001, hereafter 
IKvZOll f and lBesla fc Wul (I2007L hereafter BW07). Ob- 
served upper liraits on the g as content of debris disks 
(|Thi et al.ll200lHChen fc Jur a 2003; Robcrgc ct al. 200|) 
and theoretical inferences ([Fernandez et al.i i2006i ) sug- 
gest that gas in debris disk is non-primordial and is de- 
pleted in hydrogen. Furthermore, even whe n assum ing 
solar abundance of hydrogen, it was shown by lBWOTI that 
heating arising from H2 coUisional de-excitation, pho- 
todissociation and formation on dust, as well as cooling 
from H2 ro- vibrational lines, are negligible when com- 
pared to the photoelectric heating by grains. For these 
reasons, we do not include H2-dependent heating and 
cooling processes. We do include, however, the atomic 
hydrogen contribution to cooling via Ly a emission at 
very high temperatures (T > 5 000K). Also excluded 
from computation are coolings by vibrational/rotational 
transitions of molecular CO and CH. These molecu- 
lar species may arise from evaporating cometary bod- 
ies. However, they are believed to be quickly photo- 
dissociated aroun d the early-type st ars that are the tar- 
get of this study (|Kamp etalj 12001 . Indeed, HST UV 
observations of (3 Pic failed to detect signi ficant amounts 
of CO in the system (|Roberge et al.ll2000[ ). 

The gas temperature is solved at every spatial grid 
assuming that heating and cooling processes balance each 
other locally. 

2.1. Heating 

Relevant heating processes include: photoelectric heat- 
ing by dust (PE), photoionization of gas by the stel- 
lar light (PI), and gas-grain collisions (GG). The last 
one becomes a cooling mechanism for th e gas when 
T'gas > Tdust- For /3 Pic-like disks, IBWOTI have shown 
that photoelectric (PE) heating is by far the dominant 
heating mechanism. 

Photoelectric heating. — PE is caused by energetic stellar 
photons striking dust particles, ejecting electrons from 
the grains. Released electrons, in turn, contribute their 
kinetic energy to t he gas. The PE heating rate per unit 
volume is given bv iWeingartner fc Draing ()200 1rl as 
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Here, J'*"""' ds "^^^^"^ tts^ represents the total dust area 
per unit volume, coming from grains of radii between 



^ The IWei ngartner fc Draind 120011) formulation is a refinement 
over that in ;Drainc (19781 ). We adopt it here. 



Smin and Smax, that stellar photons can intercept, (j) is 
the charging potential of the grain, W is its work func- 
tion, Qabs is the overall absorption coefficient, Y^hi/) 
is the photoelectric yield, F^ is the stellar flux at a 
given frequency, E is energy of the ejected electrons, 
and f (E, hv) describes their energy distribution. Fol- 
lowing |^^ingartneF&|£raini ()200ll) , we take Qabs — 1- 
Throughout this paper we assume a carbonaceous com- 
position for the dust grains, giving a work function of 
l^ = 4.4eV. Also, 
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where i?high = hv — W — et/), E\ow — —ecj) and the dimen- 
sionless factor j/o is the highly uncertain photoionization 
yield. Following Weingartner fc Draine (200 li) , we take 
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silicate. (3) 



Each photoelectric electron that makes it out of the grain 
carries with it an average energy of order 1 eV. 

The charging potential (j) critically determines the heat- 
ing capability of stellar photons - only stellar photons 
with energies above a threshold {W + e0) yield PE heat- 
ing. The potential is independent of grain size and is 
obtained by equating the photoelectric charging current 
per unit area on the dust grain. 
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to the thermal electron collection current. 
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Here Sc ^ 0.5 is the electron sticking coefficient, n^ is the 
electron density, obtained in our code from the ionization 
balance computation. Tgas is the gas temperature, me 
the electron mass and fce the Stefan-Boltzman constant. 
Eq. [5] holds for positively charged grai ns (ecj) > 0), which 
is the cas e for debris disk gas den sities (IBWOTI ). Also, as 
argued bv lFernandez et all (|2006:), the charging potential 
depends logarithmically on disk parameters like electron 
density, gas temperature, etc. 

FpE depends directly on the area covered by the dust 
and through eq. [5] on electron density Uc- If we take the 
simplification that each liberated electron carries of order 
leV of energy to the gas (second integral in eq. [l]), we 
can estimate Fpo w dust area x Jo x 1 eV/e oc Uc (also see 
eq. [23]). 

Also arising from cq. [S] is the dependence of PE heat- 
ing on gas temperature, Tgas- Since e0 is generally a few 

/cbT* 3> AieTgas (BW07), e^/fcBTgas is the dominant term 

1/2 

in the equation. Hence, Jo (and in turn Fpe) oc Tgas ■ 

At higher temperatures, the Coulomb focussing is weak- 
ened and electrons are less likely to be collected by the 



positively charged grains, hence ecj) becomes more posi- 
tive, thus raising the W + ecp threshold. As a result, the 
PE heating rate is reduced at higher gas temperatures. 
There is no runaway heating. 

In practice, the double integrals of eq. [Tj are pre- 
computed for a range of values. As the integrals are 
monotonic functions of 0, the value of 4> is then obtained 
by equating the pre-computed thermal equilibrium cal- 
culation JpE to Jc. 

Gas-grain collisions. — The volumet ric heating rate due 
to GG is given by eq. [19] in lKvZOll 



Tgg = 4.0 X 10 i^ergcm-^s-i ntotaTy%^ 

X(rdust -T'gas) / TTs'^—^ds, (6) 



where ax is the thermal accommodation coefficient, 
taken to be 
cool the gas. 



taken to be 0.3. For Tgas > Tdust, gas-grain collisions 



Heating due to photoionization. — The heating rate per 
unit volume due to photoionization is 
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where Fe^^ is the energy released per unit time, per unit 
atom (i) during photoionization, the subscript j denotes 
the ionization state. The photoionization rate is precom- 
puted for a given stellar flux and recalled when needed. 
We follow atoms H through Zn and through to their sec- 
ond ionization states (see Appendix for more details). 

Heating due to Gas-Grain Drifting — Following BW07, we 
ignore drift heating due to differential velocity between 
gas and grains. This is justified in our problem where 
gas and grain are we akly co upled, and is also confirmed 
by the calculations of iGorti &: HoUcn bach (.2004) . 

2.2. Cooling 

When computing gas cooling, we include the fine struc- 
ture lines O I 63.2, 44.1, 145.5 ^m, C II 157.7^m and Si II 
34.8 ^m, the atomic cooling lines O I A6 300A and Lya, 
as well as free-free, and radiative recombination cooling. 

For gas temperatures below 5 000K, O I, C II and 
Si II fine-structure transit ions dominate the cooling 
(jHoUenbach fc McKei ll989"). Colhsionally excited levels 
decay spontaneously, converting thermal energy to pho- 
tons that are lost to space, effectively cooling the gas. 
While both collisions with free electrons and atoms/ions 
can lead to fine-structure cooling, we neglect the latter 
since electrons move faster. 

Cooling by fine- structure transitions. — Let rii be the pop- 
ulation density in level i. For sufficiently high electron 
densities (and therefore collision rates), each level is pop- 
ulated according to the local thermal equilibrium (LTE), 



g»exp(-j;i/fcBT) 
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n^, 



(8) 



where gi is the statistical weight of state i, Ei its ex- 
citation energy, nx the total number density of species 
X, and the summation runs through all energy levels. At 



lower densities, a full statistical equilibrium calculation is 
required to determine rij. The thre shold densities (ne.crit) 
between SE and LTE are listed in iHollenbach fc McKed 
(1989) for various transitions. We collect values for lines 
of interest to this study in Table [1] The typical electron 
density encountered in our systems is of order twenty or 
less. As such, we have constructed the ONTARIO code to 
perform SE calculation for all elements of interest. In 
particular, for the important cooling elements (O I, C II 
and Si II), the SE calculation is self-consistently cou- 
pled to the gas thermal equilibrium calculation. Data 
required in such calculations, as well as the numbers of 
energy levels and transitions included, are detailed in the 
Appendix. 



TABLE 1 
Relevant cooling lines 



Line 


A., 


'^c.crit 




[s-M 


[cm-3] 


CII 157.7 /im 


2.3 X 10"^ 


8.7(T/100K)"-''" 


OI 44.1 ^m 


8.9 X 10""' 


6.3 X 103(r/100K)-°'^3 


O I 63.2 ^m 


8.9 X 10-5 


6.3 X 10^{T/100K)-°°^ 


OI 145.5 Atm 


8.9 X 10""^ 


8.9 X 102 


Sill 34.8 ^m 


2.1 X 10-'' 


1.2 X 102(T/100K)°-5'' 



Note. — J^e.crlt is the critical electron density above which 
LTE is approximately correct. T is gas temperature. 

The line luminosity per unit volume for fine-structure 
de-excitation i — > j is contributed by both spontaneous 
and stimulated emission, 



^10 = [Aij 4- Bij{U^)\ Ui hv, 



tj, 



(9) 



where Aij and Bij are the usual Einstein coefficients, Ui, 
is the radiation field density at Vij. We ignore stimu- 
lated emission when computing line fluxeso However, 
we do include the stellar fluxes and dust infrared radia- 
tion when computing SE. KvZOl demonstrated that dust 
IR radiation promote population in the higher excitation 
states of O I, leading to stronger O I line fluxes. 

The energy absorbed from radiation pumping leaves 
the system in two ways. Part of it is subsequently re- 
leased as photons by spontaneous emission, part of it 
heats the gas via coUisional de-excitation. The former 
process increases the line luminosity (as calculated in eq. 
[5]), but does not contribute to the net cooling of the 
gas. As a result, we define the gas cooling rate differently 
from the cooling line luminosity (eq. [3] ) as the difference 
between coUisional excitation and de-excitation between 
levels i and j. 



Aij = {qjt n.j - q.i.j Hi) ric hv, 



I], 



(10) 



where qji and Qij are coUisional excitation and de- 
cxcitation rates, respectively. For very strong radiation 
fields in the vicinity of early-type stars, fluorescence can 
sometimes cause population inversion and lead to effec- 
tive coUisional heating of the gas. 

^ Photons arising from stimulated emission are highly 
anisotropic. In the case of stimulation by stellar photons, only 
emission from gas directly in our line-of-sight to the central star is 
detectable. Moreover, this contribution is negligible compared to 
spontaneous emission from the whole disk. We therefore decide to 
report only line fluxes from spontaneous emission. 



O I (6 300 A) and Lya transitions. — Electronic tran- 
sitions only beco me imp ortant at high temperatures 
(eqs. [34] & [35] of lKvZOlh : 
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22 800 



abundance. This is motivated below. We then com- 
pare our model output to the observed Na I profile to 
obtain values for the fitting parameters, pi sa 100 AU, 
ai = 1.0, 02 = 2.3, no = 24cm-3, a, = (p/lOOAU) x 
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Lya 



= 7.3x 10~'^ergcm 



riHI^eGXp 



— (HpAU/VSlnS. This derived distributio n is similar to 
^ /that in eq. [4] of IBrandeker et all ()2004[ ) except for two 
118 400KVljfferences: our no is a factor of 10'^ smaller as we adopt 



T 



Free-free cooling. — Free-free cooling also enters at high 
temperatures. Cooling due to e lectrons interactin g with 
ions of charge Z is (eq. [3.14] in IOsterbrock|[T989() 

Aff(Z)1.42 X lO^^Vgcm-^s-i gs Z^VTnen+, (13) 

where n+ is the ion number density and gg ~ 1.3 is the 
mean Gaunt factor for free-free emission. The total free- 
free cooling rate is the summation over all ionic species. 
The recombination cooling is similar to that in eq. [7], 
but with Fe replaced by Ae, the recombination cooling 
rate per electron for species z, and with Uij replaced by 
electron number density no: 
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Ae = EEA: 



E,ine 



3. THE B PICTORIS DEBRIS DISK 



(14) 



In order to check the reliability of ONTARIO, we ap- 
ply it to the well-studied debris disk of /3 Pictoris. The 
/3Pic system is unique in that it is the only debris sys- 
tem where both the dust and the gas components have 
been extensively mapped (see, e.g., iChen et all 120071 
and references therein). 

3.1. Input parameters 

The adopted model parameters of the /3 Pic system are 
listed in Tabled We used the same flux-c alibrated stellar 
model spectrum for /3Pic as in Fernande z et al.l f2006) , 
except we merged the model spectrum with far- and near- 
ultraviolet observations from FUSE and HST/STIS (A. 
Roberge, priv. comm.), covering 925-1 180 A and 1465- 
1 660 A, respectively, and interpolating over the gap. The 
distributio n of dust area is deduced from HST/STIS ob- 
servations (jHeap et al.l 120001 : [Fernandez et al.ll2006l ) and 
summarized by a fitting form 



0.71 



n{a )ndust 



TO exp [^{z/W) 



(15) 



where p is the cylindrical radius and z the height above 
disk midplane. The empirically determined fitting pa- 
rameters are po = 120 AU, W = 6.6(p/po)°'^^ AU, and 
To = 2 X 10"'^. We also adopt the following form to de- 
scribe the distribution of the total gas number density 
(IBrandeker et al. 2004 ): 
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(16) 

For elemental abundanc es, we assume solar abundances 
(jGrevesse &: Noel^ll993f ) except for three elements: hy- 
drogen is set to 10~^ of its solar value, helium is as- 
sumed to be zero, and carbon is set to be 20 x its solar 



ia'' hydrogen poor mix while they assumed a solar abun- 
dance of hydrogen; our vertical distribution is well ap- 
proximated by a Gaussian proflle with a full width half 
maximum (FWHM) o f H = 0.4r. This differs from 
IBrandeker et al.l (|2004D because we include the varia- 
tion of ionization fracti on with height in this work, while 
IBrandeker et al.l (|2004[ ) did not. 

The volume density of dust area (eq. [15]) vanishes 
close to the star, indicating an inner clearing of dust. 
The gas profile of the inner disk is uncertain. The Na I 
emission line can be traced in to a projected separation 
of 13 AU from the star and out to 323 AU from the star. 
Comparing the Na I seen in emission with the Na I seen 
in absorption, iBrandeker et al. ( 2004) concludes that 80- 
100 % of the circumstellar Na I seen in absorption is 
within this range. We therefore limit the gas disk to 
start at 15 AU and end at 200 AU (with more distant 
regions making negligible contributions to the total line 
luminosities). The model disk height is limited to 100 AU 
above the midplane. 

To relate the dust opacity distribution of eq. [15] to the 
dust disk mass, we assume that dust grains are spherical 
and range in sizes Smin < s < Smax, with a differential 
number distribution dn^ust/ds oc s~^'^, as characteristic 
of dust in the interstellar medium. Smax is set, some- 
what arbitrarily, to 1 cm, while Smin is set to be the min- 
imum size below which radiation pressure could remove 
the grains. For a grain initially at circular orbit, this 
occurs when the ratio of the radiation pressure to grav- 
itational force acting on the grains, /?, is equal to 0.5. 
This translates to 
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(17) 



where c is the speed of light and Qpr ~ 1 is the radiation 
pressure efficiency averaged over the stellar spectrum. 
For /3Pic, we assume iW^ = 1.75 M© and L* = IIL©. 
Taking ggiain ~ Igcm"^, we find Smin '^ 5pm. 

The infrared luminosity from the dust, integrated over 
the whole disk, is thus given by: 



-t^dust — 



L. 



At:{p'^ 



-(1 - e)7r(a^)ndust27r/9 dpdz, 



(18) 
where e ~ 0.5 is the grain reflectivity coefficient. When 
using the /3Pic dust distribution described by eq. [IB] . 

i^dust ^ 10 i*. 

To derive the grain temperature as a function of dis- 
tance from the star, we use an empirical relation derived 
from observations of the (3 Pic dust disk. From observa- 
tions (e.g. Nilsson et al. 2009) it is clear that the dust 
cannot radiate as black bodies, as that would imply too 
much sub-mm radiation. Instead, the spectral energy 
distribution of the disk is found to be well fit by a mod- 



ified black body, 



/,. 



B.{T). 



(19) 



One can attempt to relate 7 and v^ to dust properties 
(e.g. iDraind [20031 ) . but that is beyond the scope of the 
present article, where we are content that this simple em- 
pirical relation approximates the observed dust emission 
surprisingly well. Using eq. [12] for radiative equilibrium, 
where we equate the absorbed energy flux to the emitted, 
we get 
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which in turn implies 



T = 



(l-e)L.^o'c' 



32/i7r2i?2^(4 + 7)r(4 + 7) ^ 



where e is an average albedo for stellar radiation, I 
the luminosity of the star, h is Planck's constant, c 
speed of light, C is the Riemann ^-function, and . 
the (true) gamma function. Using eq. |21] to comj 
the temperatures the /3Pic dust disk, and eq. [19] 
the dust area density distribution of eq. |15| to comj 
the emission, we integrate the total emission from 
disk seen at Earth to fit 7 = 0.67 and i/q = 4 x 10^' 



by comparing to the SED summarized in iNilsson e' 
([2009) . A simplified expression for the dust tempera' 
is thus 
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0.214 
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AU 
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4 X 1014 Hz 



Returning back to the issue of elemental abun- 
dances, observations and numerical modeling of gas in 
the /3Pic system indicate that the disk may be defi- 
cient in hydrogen, consistent with a disk comprised of 
main ly metals, produced as a result of dust sublima- 
tion (iLecavelier des Etangs et al.l [20011 : iBrandeker et al.l 
12004 iFernandez et all 120061 ). In addition, the gas in 
the disk has been found to be of 20 x higher car- 
bon abundance tha n expected for solar composition 
(jRoberge et al. ' 2006) . Other observed elements are close 
to solar in relative abundance. The abundance of oxygen, 
which is important for cooling the gas, is unfortunately 
hard to constrain, beca use the available abs orption lines 
are strongly saturated (jRoberge et al.l[2006[ ). although a 
solar abundance is consistent with data. Consequently, 
for our study, we set all elemental abundances at their 
solar values except for carbon at 20 times solar, hydrogen 
at 10"'^ solar and helium abundance set to zero. 

3.2. Comparing the model to observations 
3.2.1. Disk temperature profile 

Fig. [2] shows j3 Pic disk midplane gas and dust distri- 
bution and temperature profiles together with important 
heating and cooling mechanisms. The temperature in the 
disk is determined by two major heating mechanisms: 
PE and PI heating. PE is maximum around 100 AU 
where the dust distribution peaks; PI dominates in the 
strong radiation field within ~ 25 AU of the star. With 
low gas densities gas-grain collision rates are negligible. 



TABLE 2 

/3PlC DISK MODEL INPUT PARAMETERS 



M* 


1.75 Mq 


R* 


2i?0 


L* 


11 Lq 


Disk inner cutoff 


15 AU 


Dust mass 


0.27 M® 


Dust profile 


eq. [15] 


Dust luminosity 


9.6 X 10-1 L* 


Gas mass 


8.3 X 10-1 Me 


da« nrnfilp 


Pn [Tfi] 
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Fig. 1. — Column densities for metallic gas in the /3Pic disk, 
same data as in Tab. |3] No data for O II, Na II, Ca III, Cr I are 
available as these are not observed in optical spectra. Na I and 
Fe II are strongly constrained, and no observational uncertainty 
data are available for Mg I, Al II, Si II, S I, S II, Cr II, Mn II, Fe I, 
Ni II, Zn II. 



SO that gas and dust temperatures are effectively decou- 
pled. The result is a gas temperature profile that peaks 
at the inner disk boundary and at 100 AU. The cooling is 
dominated almost exclusively by the C II 157.7 /xm fine 
structure transition. The Si II line transition becomes 
stronger in the higher temperature region around dust 
distribution peak, yet still remains '^ 100 times below 
the C II fiux. 

3.2.2. Column densities 

The edge-on orientation of the (3 Pic disk makes it pos- 
sible to observe the disk gas through absorption. In Ta- 
ble 13] we list column densities inferred from observations 
and column densities computed with ONTARIO, using the 
disk model outlined in ^ \3.1\ The same comparison is also 
presented in Fig. [1] for better visualization. 

The real gas disk is far from being cylindrically sym- 
metric, with the observed Na I profiles from the south- 
west (SW) and north-east (NE) parts of the /3 Pic disk 
differing substantially. The NE side is the brightest 
within a projected distance of 35 AU and outside 100 AU, 
and is detected to the limit of the observations at 323 AU. 
The SW side on the other hand, is slightly brighter in the 
35-100 AU region, and drops much mor e quickly outside 
100 A U than the NE side (Fig. 3 of IBrandeker et al.l 
120041) . Given that the model is based on the spatially 
resolved emission from the Na I D2 line, where the 



TABLE 3 

Column densities for metallic gas in the 13 Pic disk 



Species 


Observed 

[cm-2] 


Model 
[cm-2] 


CI 
CII 
Ol 
Oil 


(2-4; 
2 0+2-1 

(3-8; 


)x 

X 

)x 


10" 

10" 
10" 


8.9 X 10" 

6.6 X 10" 

7.4 X 10" 

0.0 


Nal 


(3.4 ±0.4) 


x 


10" 


2.5 X 10" 


Nail 








3.2 X 10" 


Mgl 


2.5 


x 


10" 


7.6 X 10" 


Mgll 
All 
AlII 
Si I 


> 2 

< 4 
4.5 

< 1 


x 
x 
x 

X 


10" 
10" 
10" 
10" 


5.3 X 10" 

1.0 X 10* 

4.5 X 10" 

1.2 X 10" 


Sin 


1 


X 


10" 


5.3 X 10" 


PI 


<7.0 


X 


10" 


2.6 X 10" 


PII 


<9.2 


X 


10" 


4.6 X 10" 


SI 


5.4 


X 


10" 


6.1 X 10" 


SII 


5 


X 


10" 


2.7 X 10" 


Cal 
Call 
Cam 
CrI 


< 2 X 10^ 
1.26tlH X 10" 


3.1 X 10* 
3.3 X 10" 

1.8 X 10" 
8.5 X 109 


Cm 


3.5 


X 


10" 


7.1 X 10" 


Mnl 


< 3 


X 


10" 


6.7 X 10^ 


Mnll 


3 


X 


10" 


4.4 X 10" 


Pel 


1 


X 


10" 


6.7 X 10" 


Pell 
Nil 


(3.7 ±0.5) 
< 7.6 


X 
X 


10" 
10" 


4.3 X 10" 
2.8 X 10" 


Nill 


1.5 


X 


10" 


2.7 X 10" 


Znl 


<7 


X 


10" 


8.9 X lOio 


Znll 


2 


X 


10" 


5.1 X 10" 



Note. — Observed column densities are compiled by 
IRoberge et al.l l[2006). Some species are unconstrained, as no ab- 
sorption line has been observed. The column density of Mg II is to 
be considered a lower limit since the line is strongly saturated. 

Na I represents only 10"'^ of the total Na, the agree- 
ment in column density between observation and pre- 
diction is rema rkable, and reinforces the conclusion by 
iBrandeker et al . (2004) that the gas seen in absorption 
is the same seen in emission. 

The model predicts most elements to be highly ionized, 
with the exceptions of O and C, due to their higher ion- 
ization potential. Only a few elements have their ioniza- 
tion fractions directly determined (as an column density 
average), but many have ionization fractions constrained 
by the data. Overall, the predicted ionization levels from 
ONTARIO are consistent with observations, with the ex- 
ception of C and S, which appear less ionized than pre- 
dicted. Interestingly, together with the elements O and 
P (which do not have their ionization level observation- 
ally constrained), these are the elements of Table [3] with 
the highest ionization potential (first ionization poten- 
tials are 10.4 [for S], 11.3 [C], 13.6 [O], and 10.5 eV [P]). 
Their ionization fractions are thus the most sensitive to 
the stellar UV spectrum at energies higher than 10 eV. A 
possibility for the discrepancy is thus that the UV flux in 
our adopted stellar spectrum is overestimated compared 
to the actual spectrum. 

There is the tantalizing possibility of doing cosmo- 
chemistry using these data. If indeed the metallic gas 
arises from grain-grain collisions, the evaporated gas 
could have the same chemical composition as the grains. 
This provides a window to study the make-up of /? Pic's 



planetesimal belt, much like what cosmochemists have 
accomplished by studying meteorites that fall on Earth. 
Among the intriguing questions to ask are: why is C so 
super-solar, and is the factor of 10 depletion observed in 
Al significant? 

3.2.3. Cooling line fluxes 

As outlined in §[131 we can estimate the total cooling 
line fluxes by estimating the PE heating rate. Assuming 
that every escaping electron heats the gas by ~leV, we 
obtain the following order-of-magnitude estimate for the 
cooling luminosity. 
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(23) 



where the factor of 2 in front of Ldust corrects for 
albedo, and R is the radius of the dust ring. This es- 
timate roughly explains the total cooling luminosity of 
~ 7 X 10"*Lq we obtain for /3Pic (Table g]), and it il- 
luminates the dependency of cooling luminosity on gas 
and dust properties. 

We compare these values to a previous model by 
IBWOTI also shown in Table H) The total cooling lumi- 
nosity there is higher by a factor of 30 and the dom- 
inant line is O I 63.2/xm line, as opposed to the C II 
157.7/.tm line here. These are explained by the three 
major differences betwee n the model we adopt here and 
that adopted by IBWOTI : 1) the g as mass (the metal- 
lic component) adopted in IBW07I is about seven times 
greater than we use here. As they also include contribu- 
tion to electron density from hydrogen, the total electron 
density is s ome fift een times greater than our value (see 
eq. [23]). 2) IBW07I used LTE to compute the level occu- 
pation, while we calculate the detailed statistical equilib- 
rium (SE). For the low electron density in the /3Pic disk 
(rio ^ 2Gcm~3 <g; ric^crit foi' some ions), an LTE treat- 
ment grossly overestimate the occupation number in the 
excited states of O I. This accounts for their strong O I 
lines. In our /3 Pic model, oxygen atoms do no contribute 
to cooling significantly. To test that there are no other 
significant differences between ONTARIO and the model 
used by B W071 . we used their gas and dust distribution 
and enforced LTE for the level populations, and con- 
firmed that these changes made ONTARIO reproduce the 
IBWOTI results. 

For the /3 Pic dust and gas profiles, turning off the dust 
IR field in our code only produces a <4 % variation on 
the O I 63.2 /im line fiux. 

Our predicted value for the C II 157.7 /im line flux 
is high enough to be detectable with heterodyne far- 
infrared spectrometer HIFI on Herschel. For a crude 
estimate of the instruments sensitivity, we use the pre- 
launch predicted sensitivity, expressed as the system 
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Fig. 2. — Physical properties of the /3Pic midplane in our modeh The upper left panel (a) shows the few important heating rates, and 
the lower left panel (c) the dominant cooling rates. The upper right panel (b) shows the gas and dust temperatures, respectively, and the 
lower right panel (d) gas and dust densities (note different scales and units for the two populations). 



TABLE 4 
( Pic cooling line luminosities 



Line 



ONTARIO output 



Luminosity 



-0J 



Flux 



[ergs 



C II 157.7 ^m 
O I 44.1 ^m 
O I 63.2 /^m 
O I 145.5 ^m 
Si II 34.8 ^m 



9.0 X 10"'* 
1.9 X 10"^^ 

7.1 X 10""' 
7.3 X 10-" 
5.6 X lO"'* 



7.2 X IQ-""' 
1.5 X 10-22 
5.7 X lO-i'^ 
5.9 X 10-18 
4.5 X lO""' 



Flux 
(no dust IR) 

[ergs 



with a similar effective integration time. 

A tentative detection of t he C II 157.7 um line was 
reported for the /3 Pic disk (jKamp et al.l 120031 ) with a 
BWfl^ of I Q-^"^ ergcm"^ s~^ {Aa result). This is more than 



cm 



7.2 X 10-15 
1.5 X 10-22 
5.5 X 10-" 
5.9 X 10-i« 
4.5 X lO-i'' 



'^'Wx stronger than our prediction and a confirmation by 
_persj;iel is necessary (we estimate it would be confirmed 
^^ wM^^'b cr confidence in ~ 6 seconds of integration time) 



X 1ft— • ■ • 

Xlgd'ortunately, the O I 63.2 //m emission predicted by 

-, CMj^ijiodel (5.7 X 10~^^ergs~^ cm"^) is far too weak to 

1 ^l^f^ -detected by PACS, which has a pre-launch estimated 
sensitivity of 10~^^ergs~^ cm"^ (5cr, Ih). The wave- 



NOTE. — Luminosities and fluxes of the important cooling lines, 
as observed on Earth, computed by ONTARIO. To estimate the 
contribution from the dust IR field, cooling line fluxes without 
dust IR field are included. 

temperature Tgys ~ 2 000K at 157.7 /im. Assuming the 
emission to be unresolved at a nearly diffraction-limited 
beam-size of 12", a flux of 7.2 x lO^^^ergs^^ cm-^ cor- 
responds to the velocity integrated antenna temperature 
0.38Kkms~^. Assuming the line to be lOkms"^ broad, 
i.e. the frequency bandwidth Av — 60 MHz, the effective 
integration time required for a 5cr-detection should be on 
the order of 



1 



6 X 107 Hz 



5 X lOkms^i X 2000 K 
G.38Kkms-i 



= 20 min. 



(24) 
The line could also be detected by PACS on Herschel, 



length region of Si II 34.8 /.tm line is not covered by any 
present or planned instrument and, as such, cannot be 
confirmed. 

4. LINE LUMINOSITIES FROM DISKS AROUND A AND F 

STARS 

Having studied the /3 Pic disk in some detail, we now 
proceed to investigate a range of potential debris disk 
configuration. We first study a fiducial disk with param- 
eters tuned to resemble the /3Pic disk of §[31 but with 
simplified gas and dust profiles. We then investigate 
the temperature and ionization rates in the disk, fine- 
structure cooling line fluxes, and identify major emission 
and absorption lines. We vary these parameters system- 
atically and report their influences on the observed cool- 
ing fluxes. In particular, we sample changes in elemental 
abundances (C, O, and H), variations in the total mass 
and distribution of gas and dust in the disk, and test a 



range of stellar spectral types. 

4.1. Fiducial disk 



TABLE 5 
Fiducial Case Input Parameters 

'mI 1.75 Mq 

-R* 2 ix0 

L* 11 L0 

Dust mass 0.5 Mgj (assuming a bulk density of 

lgcm~3 and Smax = 1cm) 
Dust profile eq. [gS]: Po.dust = 120 AU, FWHMp = 10 AU 

H/p = 0.1,'no,dust = 2.9 X lO-l*^ cm^cm-^ 
Dust size distribution dn^^i^i/ds oc s~^'^, Smin < s <lcm 
Gas mass 3.0 X 10-* Mq 

Gas profile eq. 1251 : no, gas = 24cm~'^, po.gas = lOOAU, 

FWHMp = 40 AU, H/p = 0.4 
Elemental abundances solar (Grevesse&Noels 1993), except 
[C] = 20 [C]0, [H] = 10-3 [H]q and [He] = 



value, rii — ?► gi/{go + ffi)^o- Other cooling lines come 
into importance, including the Si II 34.8 /im line, the 
Lya line, and gas-grain collisions. The O I lines, esti- 
mated to be important by BW07, are insignificant in our 
model, due to NLTE effects. 

Fig. |3] shows the gas temperature at the peak of the 
dust distribution to rise as high as 2.5 x 10'* K. Since ON- 
TARIO does not include some cooling processes becoming 
important above 5 000K (like C I A9849 C I I A2324, 
Fe II 1.26 ^m, etc. iHollenbach fc McKedri989| ). this is 
an overestimate^ There is a possibility, however, that 
the temperature is high enough for the gas to thermally 
evaporate. Let the criterion be that the sound speed 
(for Carbon) becomes comparable to the orbital escape 
velocity: 



Fiducial case input parameters are listed in table [SI 
For the gas distribution observed in the /3Pic system 
(eq. [H]), the total mass of the gas disk depends not 
only on the parameters specifying the gas profile, but 
also on the choice of inner and outer disk boundaries. 
To avoid this complication, we choose a double-Gaussian 
profile for both the gas and the dust components, 
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(25) 



Here p and z are cylindrical coordinates, Gp — 
FWHMp/v8hi2 represents the width of the radial dis- 
tribution and (Tz = (p/prcf) X FWHMz/-\/81n2 is the 
vertical profile width that scales linearly with distance. 
At Prof = 100 AU we set FWHMs, to 40 AU, giving a 
constant H/p = 0.4. The peak midplane density UQ^gas 
is set to 24 cm"^, based on the value observed for the 
/3Pic system. The interstellar medium gas density of 
6.7 X IQ-^cm-^ is chosen as a lower limit throughout 
the disk, should the Gaussian of eq. [25] fall below this 
(6.7 X 10"^ cm"^ corresponds to 1 cm^"^ if solar abun- 
dances of H and He were to be included). The Gaussian 
distribution of the dust is slightly offset from the gas in 
the radial direction (see Fig. [3]) with a steeper radial and 
vertical drop-off (FWHMp of 10 AU and H/p ^ 0.1). 

The implied dust luminosity is 1.2 x lO^'^L,^. To con- 
vert between the dust area distribution and the total dust 
mass, we computed Smin from eq. [T7] . 

4.1.1. Disk temperature profile 

Fig. [21 shows the gas and dust temperature obtained in 
our fiducial model, as well as various cooling and heating 
rates. The gas is thermally decoupled from the dust, as 
the low gas-grain collision rate is not able to equilibrate 
the two population. We find that PE heating dominates 
the gas heating where dust density is high, and that PI 
heating dominates where dust density is relatively low. 
Over most of the disk, cooling is dominated by the C II 
157.7 /xm line. However, as gas temperatures rises above 
~ 100 K, the cooling fiux of C II saturates as the up- 
per occupation number ni approaches the LTE limiting 
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(26) 



As this temperature limit is below our (over-)estimate, 
we cannot rule out gas evaporation from the disk; a more 
detailed study, outside the scope of the present paper, 
would be required. 

4.1.2. Cooling line fluxes 

Table [Bj lists cooling line luminosities integrated over 
the whole disk for the fiducial case. C II 157.7 /xm re- 
mains the most luminous line, with Si II 34.8/im follow- 
ing within an order of magnitude, while O I lines are very 
unde rluminous compared to previous study (|Besla fc Wul 
I2007D due to NLTE effects. Most of the line flux is gener- 
ated in region of highest gas and dust densities: 93 % of 
the total C II 157.7 /^m flux arises in the region extending 
radially from 100 to 140 AU and vertically up to 20 AU. 

For C II, which is close to LTE, the line flux is domi- 
nated by coUisional processes (eq. [S] ~ eq. [TO]), which 
in turn is determined by the heating processes (in this 
case, photoelectric heating). Compared to the /3 Pic case, 
the fiducial case has a higher gas temperature in the 
dust maximum and consequently a smaller photoelectric 
heating rate. This explains largely the lower C II flux 
in the latter. This difference is further compounded by 
enhanced photoionization heating in the /3 Pic case due 
to a larger amount of gas closer to the star. O I, on 
the other hand, experiences an electron density that falls 
far below its critical density for LTE and its line flux 
arises largely from excitation by stellar photons (radia- 
tive pumping) and not collisional excitation which leads 
to cooling. This is reflected in the third column of Ta- 
ble [6] which shows that the cooling fraction (fraction of 
the line flux responsible for cooling the gas, the ratio of 
eq. [TO] to eq. [H]) for O I is 0.06 < 1, while it is 0.96 for 
GIL 

4.1.3. Ionization profile 

In our fiducial model, the metallic gas is strongly ion- 
ized by the central star. Figs. [3]&[n] show the neutral 
fractions of Ga, Na, and G, the main contributers to the 
electron density. We also show the fractions of Ga II and 
Ga III since a significant fraction of Ga is ionized a second 
time. 

•^ However, our tests indicate that at the low gas density assumed 
here, these cooling mechanisms are negligible even at temperatures 
as high as 10^ K. 
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Fig. 3. — Similar to Fig. [2] but for the fiducial model . Due to the higher "dust/^gas value at the dust peak, the gas temperature there 
experiences a spike and rises beyond 5 000 K (see ii l4.1.1l for discussion) . 
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TABLE 6 

Fiducial case cooling line luminosities 

Luminosity Flux at 20 pc Cooling 



[L 



01 



[erg s -"- 1 



Fraction^ 



0.96 

1 
0.06*= 

J~ 
0.63 



L CJ - - 

ClI 157.7 /im 2.3 x 10"** 1.8 x 10"^'' 

Ol44.1/.4m 3.5x10-1'' 2.8 x 10"^^ 

Ol63.2/.4m 1.5x10-1° 1.2x10-" 

Oll45.5Mm 1.4x10-" 1.1 x lO-^* 

SiII34.8/jm 3.2x10-^ 2.6 x 10-^^ 

Note. — Cooling line luminosities and fluxes measured at 20 pc, 
as predicted by ONTARIO for the fiducial model. 

^Fraction of line luminosity from coUisional processes (the ratio 
of eq. 1101 to eq. [9]); the remainder is due to radiative pumping 
and does not contribute to cooling of the gas. 

''Total O I cooling fraction summed over all 3 transitions. 

The ionization profile is explained by a competition be- 
tween stellar pliotoionization and recombination. Rate 
for the former falls off as one moves away from the star, 
while rate for the latter scales with electron density (and 
weakly depends on temperature through the recombina- 
tion coefficient). So the neutral fraction is the lowest 
near the star and it rises upward until it reaches a maxi- 
mum around 100 AU. Outward of this distance, the sharp 
drop of gas density reduces the recombination rate faster 
than the decrease of the ionization rate, resulting again 
in a more ionized gas. The neutral fractions exhibit a 
narrow dip at around 120 AU. This is due to the temper- 
ature peak (as shown in Fig. [S]), which in turn is due to 
the dust peak at that location. As is seen in Fig. [S] the 
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Fig. 4. — Midplane neutral (unbroken line) and ionization 
(dashed for one time ionized, and dash-dotted for two times ion- 
ized) fractions of three selected elements: Ca, Na, and C. 



neutral fractions decrease monotonically away from the 
midplane due to the reduction in recombination in lower 
densities. 



11 




20 40 60 

Vertical height z [AU] 



100 



Fig. 5. — As in Fig.|4l but for the vertical distribution at 100 AU 
(thick line) and 120 AU (thin Une). 



4.1.4. Optical/UV absorption and emission lines 

Observationally, an advantage of the IR cooling lines 
is that the stellar photosphere is relatively dark at those 
wavelengths, making faint emission lines more easily de- 
tectable. Should the disk be spatially resolvable or be 
observed edge on, strong transitions in the optical/UV 
might still be easier to detect. Inde ed, gas around 
;gPic was first de tected in absorption (jSlettebakl 119751 : 
iHobbs et al.lll985l ). and then in spatially resolved emis- 
sion f rom light scattered in atomic resonance lines (Na I 
D2,i; lOlofsson et al.ll200ll ): apart fr om the tentative de - 
tection of C II 157.7 fj,m. emission (jKamp et al.l 120031 ) . 
cooling lines have yet to be observed from the disk. With 
this in mind, we have computed equivalent widths for ab- 
sorption lines assuming an edge-on disk (Table Uil , and 
total luminosities for light scattered in lines with strong 
transitions (Table [H]) . Resolved observations of the Ca II 
K absorption line around /3Pic (jCrawford et al.llI998D 
show it to be ^ 2kms~^ wide. We assume absorption 
from the 'stable component' of other species to be sim- 
ilarly broadened, resulting in increasing saturation for 
lines with equivalent widths approaching 2kms~^. 

4.2. Parametric survey 

To detect circumstellar gas in absorption requires a 
special geometry, and to detect light scattered in gas 
most likely requires a resolved disk. A comprehensive 
survey for tenuous gas in debris disks is therefore best 
served by the IR cooling lines. In this section, we study 
how the cooling line luminosities depend on various disk 
parameters, focusing on the C II 157.7 /im and O I 
63.2 /im lines, since these fall in the wavelength window 
of the Herschel telescope. The parameter study is ob- 
tained by modifying one parameter at a time while keep- 
ing all other parameters at their fiducial values. Disk 
parameters and their ranges (provided for easy reference 
in Table H]) that we cover include: 



TABLE 7 
[JCiAL Case Dominant Absorption Lines 


Species 


'^vacuurr^ Eq 

fAl 


Width'' 
[mA] " 


NaT 


5891.58 


0.86 


NaT 


5897.56 


0.43 


Call 


3934.78 


OT 


Call 


3969.59 


OT 


Till 


3350.37 


0.82 


Till 


3362.18 


0.59 


Till 


3373.77 


0.51 


Till 


3384.74 


0.39 


Mnll 


3442.97 


0.15 


Pel 


3441.59 


0.31 


Fel 


3720.99 


0.88 


Fel 


3735.93 


0.36 


Pel 


3861.01 


0.53 



''Only lines with A > 3300 A are presented. 

''OT - optically thick line, defined to have an equivalent width 
that corresponds to a velocity > 2knis~-'. 

1. peak gas density. Range covered is from 0.01 /3Pic 
density (2.4 x 10"^ cm'^) to -- 2.7 x 10^ cm~3, the 
latter corresponds to a column density above which 
the C II 157.7 /im line becomes radially optically 
thick. For the spatial distribution we adopt, these 
densities correspond (through eq. [5S] with fiducial 
parameters) to a total gas masses of 2.1 x lO^^Af0 
and 3.3 X 1O~^M0, respectively. 

2. total dust mass. We sample dust mass from 1 zo- 
diac (taken as M^ust = 10~^ M^; with fiducial disk 
parameters, this corresponds to Ljust/^* ^ 10~^) 
to 3 Mq , at which point the dust becomes radially 
optically thick (with Ldust/-^* ^ 3 x 10~'^). 

3. spatial distributions of the gas and dust compo- 
nents, including location of peak density, radial 
FWHM and vertical scale height, spanning about 
an order of magnitude in each parameter. 

4. elemental abundances. We vary carbon and oxygen 
abundances, respectively, from solar to 20 x solar, 
to account for the type of abundance anomaly de- 
tected in the (3 Pic disk. 

5. stellar spectral type. We study main sequence stars 
with effective photosphere temperatures ranging 
from 6 500K to 10 000 K. For cooler stars, molec- 
ular chemistry is likely to be important, rendering 
the ONTARIO code invalid. We interpolate between 
the ZAMS models i n Table 15 .14 in Allen's Astro- 
physical Quantities ()Coxll2000l ) to obtain the stellar 
luminosities and sizes. 



4.2.1. Gas and dust densities 

We show the integrated luminosities from the C II 
157.7/im and the O I 63.2 /im lines as functions of gas 
peak density and dust mass in Figs. [6] & [7l respec- 
tively. In systems with little dust, gas is primarily heated 
through photoionization, so the C II luminosity scales 
almost linearly with the amount of gas. As more dust 
is added, photoelectric heating takes over and the C II 
flux starts to increase with both dust and gas density 
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TABLE 8 
Fiducial Case Dominant Emission Lines 



Species 


^vacuun^ 


Luminosity 


Flux at 20 pc 








Lq 


[ergs- 


1 cm-2] 


CI 


9826.8 A 


1.5 


X 10"* 


1.2 X 


10-1^ 


CI 


9853.0 A 


4.4 


X 10"* 


3.5 X 


10-15 


CII 


157.7^111 


2.3 


X 10-8 


1.8 X 


10-15 


Ol 


63.2/jm 


1.5 


X 10-10 


1.2 X 


10-17 


Nal 


5897.6 A 


7.6 


X 10-8 


6.1 X 


10-15 


Nal 


5891.6 A 


1.5 


X 10-^ 


1.2 X 


10-1-1 


All 


3945.1 A 


2.2 


X 10-10 


1.8 X 


10-17 


All 


3962.6 A 


4.4 


X 10-10 


3.5 X 


10-17 


Sill 


34.8 /im 


3.2 


X 10-9 


2.6 X 


10-16 


SI 


25.2^m 


2.4 


X 10-10 


1.9 X 


10-17 


Call * 


3934.8 A 


9.4 


X 10-'^ 


7.6 X 


10-13 


Call * 


3969.6 A 


7.0 


X 10-'^ 


5.6 X 


10-13 


Call 


7293.5 A 


3.4 


X 10-^ 


2.8 X 


10-1-1 


Call 


7325.9 A 


3.2 


X 10-^ 


2.5 X 


10-14 


Till 


4534.5 A 


4.2 


X 10-10 


3.4 X 


10-17 


Till 


4983.1 A 


4.6 


X 10-10 


3.7 X 


10-17 


Till 


3350.4 A 


1.6 


X 10-^ 


1.3 X 


10-14 


Till 


3362.2 A 


1.1 


X 10-^ 


8.8 X 


10-15 


Till 


3373.8 A 


1.1 


X 10-^ 


8.8 X 


10-15 


Till 


3384.7 A 


8.4 


X 10-8 


6.8 X 


10-15 


Cm 


8002.3 A 


7.4 


X 10-8 


5.9 X 


10-15 


Cm 


8127.5 A 


5.8 


X 10-8 


4.7 X 


10-15 


Cm 


8231.9 A 


4.7 


X 10-8 


3.8 X 


10-15 


Mnl 


4034.2 A 


2.2 


X 10-9 


1.8 X 


10-16 


Mnl 


4031.9 A 


3.0 


X 10-9 


2.4 X 


10-16 


Mnll 


3443.0 A 


3.5 


X 10-8 


2.8 X 


10-15 


Mnll 


3461.3 A 


1.8 


X 10-8 


1.5 X 


10-15 


Pel 


3721. OA 


1.4 


X 10-^ 


1.1 X 


10-1-1 


Pel 


3735.9 A 


8.6 


X 10-8 


6.9 X 


10-1-1 


Pel 


3821.5 A 


9.9 


X 10-8 


7.9 X 


10-15 


Pel 


3861. OA 


1.3 


X lO-'' 


1.1 X 


10-1-1 


Pc II 


12570.9 A 


1.6 


X 10-^ 


1.3 X 


10-1-1 


Pell 


13209.9 A 


4.3 


X 10-8 


3.5 X 


10-15 


Pell 


16440.0 A 


4.1 


X 10-8 


3.3 X 


10-15 


Pell 


25.99^*111 


6.1 


X 10-8 


4.9 X 


10-15 


Nill 


6668.6 A 


1.9 


X 10-8 


1.5 X 


10-15 


Nill 


7379.9 A 


8.0 


X 10-8 


6.4 X 


10-15 


Nill 


7413.7 A 


3.1 


X 10-8 


2.5 X 


10-15 



Note. — The two Ca II lines indicated by {*) are optically thick. 
While our code includes the scattering of stellar flux from the star 
to the computational bin, it does not account for the scattering of 
gas emissions along the line of sight of the observer. Therefore, the 
computed values may significantly overestimate the luminosities in 
the 2 saturated calcium lines. 
''Only lines with A > 3300 A are presented. 
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Fig. 6. — C II 157.7 ^m integrated line luminosity as a function 
of peak gas density and total dust mass. Contours are separated 
by a factor of 10. The arrow shows the direction of increasing 
temperature. 




10 10 

Gas peak density [cm"' 

Fig. 7. — O I 63.2 ^m integrated line luminosity as a function of 
peak gas density and total dust mass. A 10-level O I atom is used 
in the computation . The line luminosity is largely independent of 
dust mass and scales roughly linearly with gas density. 



TABLE 9 

Parametric survey: overview of sampled parameters 



Parameter sampled 




Min 


Max 


Gas peak density 


2.4 


X 10-1 cm-3 


2.7 X 103cm- 


(total gas mass) 


(2.1 


X 10-6 A%) 


(3.3 X 10-2 A/( 


Dust mass 




10-3 Me 


3Mg5 


Gas peak location 




30 AU 


190 AU 


Gas FWHM 




5AU 


95 AU 


Gas scale height 




0.05 


0.95 


Dust peak location 




30 AU 


190 AU 


Dust FWHM 




5AU 


95 AU 


Dust scale height 




0.05 


0.95 


Spectral types (T^a) 




6 500K 


10 000 K 


Elemental abundances 




solar 


20 X solar 



again loses its dependence on the dust mass. Despite mi- 
=SB£ features, one can summarize the results in Fig. [5] as 



-^-21i2Eat the C II flux rises roughly linearly with gas density. 



(eq. [23]). For very massive dust disks, gas temperature 
is so high that the C II flux becomes saturated and once 



Except for the most massive dust disks where SE be- 

„ gins to be dominated by the strong dust IR field, the con- 

33 tours for the O I 63.2 /xm line flux appear independent 

]^9 of the dust mass for most of our disks (Fig. [T]). This is 

19 explained by the fact that the line flux is almost com- 

33 pletely dominated by radiative pumping, which masks 

19 any features from the coUisional processes. To remove 

19 this masking we repeat the computation but with the O I 

8 atom having only the three lowest energy states (as op- 

pos ed to the 10 levels in our normal calculation, ). Data 

both with and without the dust IR field are presented 

in Figs. [8] andini respectively. With radiative pumping 

disabled, the results show similar dust mass dependence 
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Fig. 8. — Similar to Fig. [7]except that here a 3-level O I atom is 
adopted. This removes the possibiUty of radiative pumping. Line 
luminosities are much weaker and they depend on dust mass. 



IR field due to varying dust profiles has no effect on the 
cooling line fluxes. 

4.2.2. Elemental abundances 

Since carbon is the dominant electron donor in the 
fiducial disk, increasing the carbon abundance increases 
He- This, in turn, enhances the photoelectric heating as 
is demonstrated by eq. [5] . The effects on the line fluxes 
are seen in Fig. [TT] (and Table [TU|. In contrast, increasing 
the O abundance has no significant effect on the thermal 
balance, besides from raising the O I 63.2 /im fiux by 
providing more atoms to be radiatively pumped. 

We assume the gas to be poor in hydrogen. Raising 
the hydrogen abundance to its solar value increases the 
gas-grain collisions which cools the gas. This lowers the 
equilibrium gas temperature and slightly increases the 
photoelectric heating (and the corresponding C II cool- 
ing flux). But the effect is not significant. In fact, one 
can remove all elements from the gas except for C, O, 
Si without causing much change in the cooling line lumi- 
nosities (Table [TT 




TABLE 10 
Line luminosities as functions of elemental abundances 



10 10 10 

Gas peak density [cm~ 



Fig. 9. — Similar to Fig. [7]except that here a 3-level O I atom is 
adopted, and dust IR field is disabled. 



as that of the C II line. For more massive disks the flux 
from the cooling lines is dominated by the population 
equilibrium determined by the dust IR field. 

In Fig. llOi cooling luminosities are presented as a func- 
tion of a number of parameters that define the geometry 
of the gas and dust distributions. 

Overall, the C II 157.7 /im luminosity is at the greatest 
when peaks of the gas and dust distributions coincide, al- 
lowing for maximum photoelectric heating, or when bulk 
of the gas is located close to the star, allowing for max- 
imum radiative pumping and photoionization heating. 
The O I luminosity steadily increases as the gas distribu- 
tion is moved closer to the star as this line is dominated 
by radiative pumping. In contrast, the radial width and 
the vertical scale height for the distributions do not affect 
line luminosities appreciably. Besides having a minimal 
upregulating effect on the O I 63.2 /im line fiux (as ob- 
served for the /3Pic configuration), changes in the dust 



Modification'' 



Fiducial 
"C = "0,0 
no = 20 no,0 
"H = "H,0 
C, O, Si only^ 



C II 157.7/im 

[Lq] 



2.3 X 10"* 
3.8 X 10-1° 
2.3 X 10"* 
2.6 X 10-8 
2.3 X 10-*^ 



O I 63.2 /im 

\Lq] 



1.3 X IQ-io 

1.2 X 10-1° 
2.5 X 10-^ 

1.3 X 10-10 
1.3 X 10-1° 



^Fiducial case inputs are used except for the parameter men- 
tioned. 
I'All elements except for C, O, Si and H are removed. 



4.2.3. Stellar spectral type 

In Fig. HH C II 157.7 Aim, O I 63.2 /xm and 145.5 /im 
luminosities are plotted against stellar spectral types 
(represented by their respective effective temperatures). 
The C II 157.7 /im line remains dominant for all spec- 
tral types. Its luminosity is constant for T^g < 8 000K. 
Around more luminous stars the excited level of C II is 
populated mostly by stellar radiation (through fluores- 
cence of UV lines), rather than by collisions with elec- 
trons. In this regime, the C II 157.7 /im flux is greater 
than the actual cooling and it increases with TcS until 
Toff reaches ~ 10 000 K. Relatedly, level population in 
the excited states of the O I atoms is dominated by ra- 
diative pumping for all stars with Tcff > 7 500K. As a 
result, the line flux rises with Toff above 7 500K. 

5. CONCLUSIONS 

We have produced a thermal and statistical equilib- 
rium code to speciflcally address gas emission from dusty 
gas-poor disks around A-F stars, motivated by the recent 
launch of the far-infrared observatory Herschel, as well as 
the discovery of gas in a number of debris disk systems. 
We study the effects of, among other things, photoelec- 
tric heating from dust grains, photoionization heating of 
the gas, thermal equilibrium reached in the disk, statis- 
tical equilibrium of atoms in the disk, and cooling by 
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Fig. 10. — C II 157.7/im and O I 63.2 fim integrated line luminosities as functions of gas (top) and dust (bottom) peak cylindrical radii 
(left panels), their radial FWHM (middle panels) and their scale heights (right panels). One parameter is varied at a time. Only the left 
panels show strong dependence - the line radiation is enhanced whenever the peak radii for the dust and the gas coincide, or when bulk of 
the gas lies very close to the star. 
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Fig. 11. — C II 157.7/im and O I 63.2 fj,m integrated line luminosities as function of C (left) and O (right) abundances. The former rises 
with carbon abundance, while the latter with oxygen abundance, but the reasons behind are somewhat different (see text). 



infrared atomic lines. Details about the code (Ontario) 
are presented in the appendix. 

Using this code and disk parameters that resemble the 
observed /3Pic gas disk, we have computed a range of 
observables, including the line luminosities from infrared 
transitions, the emission and absorption of optical and 
UV lines, and the column densities of metals in the disk. 

As a guidance for observations to be carried out on 
Herschel and other telescopes, we explore the dependence 
of the infrared line luminosities on disk parameters. Our 
findings are summarized below: 



• for most of the configurations surveyed, we find 
that the C II 157.7 /im line is the most luminous. 
The expected flux for the /3Pic disk is ^ 10^^^ 
ergs~^ cm~^, which should be easily detectable by 
Herschel. 

• the O I lines, initially thought to be important cool- 
ing lines (Besla_& Wu 2007), have underwhelming 
luminosities due to NLTE effects. Most of the lu- 
minosities in these lines arise from fluorescence of 
stellar UV photons, as opposed to from collisional 
cooling of the gas. The flux from the O I 63.2 /itm 
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Fig. 12.— C II 157.7 Atm and O I 63.2/im,145.5Atm integrated 
line luminosities as function of effective temperature (T^ff) of the 
central main-sequence star. Around stars with T^ff > 8 500K, the 
level population of C II is dominated by radiative pumping, while 
population in O I is dominated by radiative pumping for stars with 
T^fi > 7 500K. 



line will be too faint to reasonably be detected by 
Herschel 

• over the parameter range we sample, the C II line 
flux scales roughly linearly with the gas density, 
with the carbon elemental abundance, and rises 
with the effective temperature of the star. Line 



fluxes also reach maximum when the dust and gas 
components in the disk align radially. 

We note two major caveats in our model: 

• in some of our models, the gas temperature rise 
above 5 000K. This is possibly an artificial result 
because we have not included all relevant gas cool- 
ing mechanisms. We have performed preliminary 
tests by consi dering other fine transition lines listed 
in (|Hollcnbac h fc McKed HQSQ) but do not find 
them of importance. However, if this high tem- 
perature feature does occur in real disks, it implies 
loss of metallic gas by thermal evaporation. 

• we have assumed that molecular species are unim- 
portant for the thermal equilibrium of the gas. For 
cooler stars, a significant amount of the gas may 
be in molecular phases (CH, CO,...) and this may 
impact our predictions of infrared line luminosities. 

By comparing Ontario predictions with observations, 
one may hope to infer physical properties of the de- 
bris disks, including gas density and elemental abun- 
dances. With the new observing windows being opened 
up by Herschel, SOFIA, and ALMA, we expect many 
new observational results that will require interpretation 
by models such as Ontario. 
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the manuscript. 



APPENDIX 

THE ONTARIO CODE 

An overview of Ontario's program fiow is shown in Fig.[T31 Ontario in principle consists of a master routine and 
four principal subroutines. The four subroutines are: 

1. Geometry solver: defines the geometry and produces observables from the solution. Our models are currently 
two dimensional (radius and height), but could in principal be of any geometry. The geometry solver computes 
the projected flux from the defined geometry. 

2. Ionization balance: given a temperature, radiation field, and elemental abundances, it guesses an electron density 
and then iterates until the derived electron density (computed from the ionization of the actual elements input) 
agrees with the assume d electron den sity. Ionization cross-sections and recombination coefficients are retrieved 
from the CLOUDY code ()Ferland ll2005'l. and are restricted to the first 30 elements (atomic numbers 1-30, hydrogen 
through zinc). To estimate the radiation field, we use NextGen stellar atmospheres (Hauschildt ct al. 1999) with 
solar abundance and main sequ ence surface logq. No chrom osphere emission is estimated in general, but the 
inter stellar radiati on field f from IWeingartner fc Drainell200lT ) and cosmic ray ionization (2 x 10~^^s~^ atorn"^, 
from lSpitzerill978f ) are included. 

3. Thermal balance: given an ionization state, dust density and radiation field, this routine computes all relevant 
heating and cooling mechanisms (as outlined in §[2|) assuming an initial temperature, and then modifying the 
assumed temperature until the heating equals the cooling. If the resulting temperature is different from what 
was assumed for the ionization balance, that routine recomputes the ionization state, which is then input to the 
thermal balance routine for a new thermal balance computation. 

4. Level population: given a radiation field, an electron density and electron temperature, the electronic energy level 
population is calculated using a statistical equilibrium (SE) of all considered transitions for given elements. This 
routine is mostly used at the end, when the physical state of the disk is already computed, as most levels and 
species do not participate in the cooling of the disk. The exception is O I, C II, and Si II where we do include a 
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Fig. 13. — Model flowchart describing the ONTARIO code. The dashed path iUustratcs how the code is run when fitting observations. 



statistical equilibrium computation already in the thermal balance computation, because of the important cooling 
lines of those spec ies. The radiative data for all atomic species was t aken from the NIST Atomic Spectra Database 
()Ralchenkoll2005| ) and the Kurucz Atomic Spectral Line Database (|Kurucdil995 ). The transition data from the 
two sources were combined to produce a more complet e energy level s et. C ollisional coefficients f or the cooling 
C II and O I lines were obtaine d from the Iro n Projec t (lPradhanll2000( ) and lSilva fc Viega^ (|2001l ). respectively; 
for C I, Si I, Si II and S I from lHollenbach"fc McKee (|1989f) . Collisional data for other electron transition lines 
were taken from the Iron Project ( PradhanI 12000 ) for C II, Ni II, Fe I, Fe II. At higher temperatures, where 
there were no data available, the collisional strength was extrapolated using a first order polynomial fit in the 
log-log space. The number of levels and transitions for different elements included in Ontario are summarized 
in Table [nj 

ONTARIO assumes the disk is optically thin in the continuum, so that simplified radiative transfer can be used for both 
the ionization and thermal balance, which simplifies the computation greatly and allows for more general geometries 
(compared to, e.g., cloudy). 

Collisional strengths 7 for the electronic transition lines for the majority of atomic species are not available from 
literat ure. 7 is calculated by taking into account a number of quantum effects (see, for example, I Zhang fc PradhanI 
(|1995f )) and no simple relationship exists between the collisional strengths and the energy associated with a transition 
line. Since we were only able to obtain collisional data for the Si II fine structure cooling line and not other higher level 
electronic transitions, we investigated whether this might result in incorrect determination of the gas temperature. 
Results of running the code with and without C II non-cooling collisional strengths (that is, setting qij = qji = for 
i, J 7^ 1, 0) for a number of bins produced no noticeable change in the computed temperatures (within 1 K). Therefore, 
all non-cooling electron transitions can be treated as purely radiative for gas densities representative of debris disks. 

SE calculation is the most computationally expensive part of the numerical model. An example on the extreme side 
is Fe II, which has 790 electron energy levels and 54054 transitions. To solve SE for this atomic species matrix of 
54 054 elements has to be inverted for each bin in the grid (the actual number of levels used was reduced to 300 lower 
levels: the population of the higher levels is negligible, and this reduction significantly improved computation time). 
Furthermore, while O I, C II and Si II have substantially less electron energy levels (150, 191 and 143 respectively) 
and corresponding transitions (1048, 1347 and 927 respectively), SE for these species participating in cooling of the 
gas has to be solved multiple times for each bin when recursively computing thermal equilibrium. 

To avoid a full SE computation, 2 and 3-level atom a pproximation for C II and O I fine structure line cooling , 
respectively, were previously used in similar simulations (iHollenbach fc McKed [19891 iKamp fc van Zadelhofil [20011) . 
We used Ontario to test the validity of this assumption for the hotter stars sampled by our model. For C II around a 
Tcff = 10 000 K star we found a factor of 4 increase in the cooling line flux when using an ion with all the energy levels 
compared to a 2-level approximation. The difference was determined to be due to fluorescence (radiative pumping to 
higher levels). While the — >■ 1 energy transition is radiatively forbidden, ground state electrons can get photo-excited 
to higher energy levels and then spontaneously decay to the first excited level. Our results show that a 2-level C II 
approximation underestimates the population of level 1, and consequently the 1 — >■ emission fiux. The effect is more 
pronounced closer to the star where the photo-excitation is more important. Similar results were seen for O I and Si II. 
Furthermore, since gas densities are substantially below no,crit for these species, the effect of radiative pumping is much 
stronger. In fact, using a significantly reduced number of levels for O I and Si II can for some system configurations 
underestimate flux by more than a factor of a 100. In addition, strong radiative pumping can overpopulate the flrst 
excited energy state beyond its LTE level, leading to stronger collisional de-excitation than excitation flux. In this 
regime the flne-structure transition becomes a heating line, injecting stellar energy into the gas. This effect is absent 
when radiative transitions are excluded. 
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We further tested the effect of reducing the number of energy levels for the cooling atoms and determined that a 
limited reduction did not significantly effect the computed temperature and line fluxes. Using 10 level approximations 
for C II, Si II and O I species kept temperature close to the full-level species results and produced similar line fluxes. 
However, the decrease in computation time was signiflcant. For example, reducing C II from 191 to 10 energy levels 
decreases execution time by a factor of ^ 365 (the square of the level reduction). Therefore, we used these limited 
reduction cooling species approximations in all further numerical runs. 
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TABLE 11 

Number of energy levels and transitions included in Ontario 

Species # of Energy Levels # of Transitions Refs for coUisional rates. 

Hollcnbach fc McKee (1989) 
Pradhan (2000) 
Silva fc Viegafi (2001) 



HoUenbach & McKec (1989) 
Hollenbacli & McKec (1989) 
HoUenbach fc McKee (1989) 



Pradhan 
Pradhan 
Pradhan 



Note. — The radiative data for all atomic species was taken from the NIST Atomic Spectra Database IIRalchenkoll2005l) and the Kurucz 
Atomic Spectral Line Database IIKuruczlll995l) . No coUisional data were available for Na I, Mg II, Al I, Al II, Ca II, Ti I, Ti II, Cr II, Mn I, 
Mn II, Zn I and Zn II. 

Note. — 
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